#!/usr/bin/env python
# vim: set fileencoding=utf-8 :
#
# This file is part of Coruja-scripts
#
# Coruja-scripts is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 3 of the License.
#
# Coruja-scripts is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with Coruja-scripts.  If not, see <http://www.gnu.org/licenses/>.
#
# Copyright 2011 Grupo Falabrasil - http://www.laps.ufpa.br/falabrasil
#
# Author 2011: Pedro Batista - pedosb@gmail.com

import argparse
import sys
import numpy as np
from config import *
from HTKParameters import HTKParameters as hpm

class MatchException(Exception):
	pass

class Dtw():
	class Dict(dict):
		def __getitem__(self, y):
			try:
				return super(Dtw.Dict, self).__getitem__(y)
			except KeyError:
				return np.inf

	def __init__(self, s1, s2, r=None):
		self.s1 = s1
		self.s2 = s2
		self.r = r
		self.M = 0

	def d(self, i, j):
		fi = self.s1.get_frames(i)[0]
		fj = self.s2.get_frames(j)[0]
#		print fi, fj
		return sum([(fi[k] - fj[k]) ** 2 for k in xrange(len(fi))])

	def lp(self, i, j):
		dij = self.d(i, j)
		return min([
			self.D[i-1, j] + dij,
			self.D[i-1, j-1] + 2 * dij,
			self.D[i, j-1] + dij])

	def p(self, i, j):
		dij = self.d(i, j)
		p1 = self.D[i-1, j] + dij
		p2 = self.D[i-1, j-1] + 2 * dij
		p3 = self.D[i, j-1] + dij
		#Rabinner Eq. 4.154
		if p1 > p2 or p3 > p2:
			self.M += 1
		else:
			self.M += 2
		return min(p1, p2, p3)

	def match(self):
		#sakoeChiba78 Fig.4
		i = 0
		j = 0
		self.D = self.Dict({(0,0) : 2 * self.d(0, 0)})
		while True:
			i = i + 1
			if self.r is not None:
				if i >= self.s1.n_samples:
					jr = int(j - self.r)
					i = 0 if jr < 0 else (jr if jr < self.s1.n_samples else None)
					j += 1
					if j >= self.s2.n_samples:
						break
			else:
				if i >= self.s1.n_samples:
					i = 0
					j += 1
					if j >= self.s2.n_samples:
						break
			if i is None:
				raise MatchException('Error: cannot reatch the final state\n\t Incriase r?')
			self.D[i, j] = self.p(i, j)
#			print i,j,self.D[i,j]
		return self.D[self.s1.n_samples-1, 
				self.s2.n_samples-1] / self.M #self.s2.n_samples

def dtw_distance(s1, s2):
	if not isinstance(s1, hpm):
		s1 = hpm(s1)
	if not isinstance(s2, hpm):
		s2 = hpm(s2)
	return Dtw(s1, s2, r=config.getfloat('dtw', 'r')).match()

if __name__=='__main__':
	parser = argparse.ArgumentParser(description="Compute accumulated DTW distance.")
	parser.add_argument('-1', '--signal-1', required=True,
			help='First signal')
	parser.add_argument('-2', '--signal-2', required=True,
			help='Second signal')
	parser.add_argument('-r', '--window', type=int, default=None,
			help='Adjustment window (Sakoe e Chiba 78 - III.B)')
	args = parser.parse_args()
#	print Dtw(hpm(args.signal_1), hpm(args.signal_2), args.window).match()
	print dtw_distance(args.signal_1, args.signal_2)
